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ABSTRACT 


Recent chemical abundance measurements of damped Lya absorbers (DLAs) re¬ 
vealed an intrinsic scatter in their metallicity of ^ 0.5 dex out to z ~ 5. In order 
to explore the origin of this scatter, we build a semi-analytic model which traces the 
chemical evolution of the interstellar matter in small regions of the Universe with 
different mean density, from over- to underdense regions. We show that the different 
histories of structure formation in these regions, namely halo abundance, mass and 
stellar content, is reflected in the chemical properties of the protogalaxies, and in par¬ 
ticular of DLAs. We calculate mean metallicity-redshift relations and show that the 
metallicity dispersion arising from this environmental effect amounts to ~ 0.25 dex 
and is an important contributor to the observed overall intrinsic scatter. 

Key -words: ISM: abundances, evolution, Galaxies: star formation, quasars: absorp¬ 
tion lines, abundances 


1 INTRODUCTION 

Damped Lya absorbers (DLA), defined as quasar absorp¬ 
tion systems with neutral column densities Nhi > 2 x 
cm“^ (Wolfe et al. 1986), dominate the neutral gas con¬ 
tent of the Universe in the redshift range z = 0 — 5 
and are likely the progenitors of low redshift galaxies 
(Wolfe, Gawiser & Prochaska 2005). The chemical proper¬ 
ties of DLAs can be determined with great precision, and 
provide a unique probe of the properties of cold neutral gas 
out of which stars form at high redshifts. 

Numerous surveys conducted over the past several years 
have provided extensive samples of high-redshift absorbers 
(e.g. Peroux et al. 2005; Prochaska, Herbert-Fort & Wolfe 
2005; Noterdaeme et al. 2012). High resolution observations 
of nearly 250 of these systems established a statistically sig¬ 
nificant decrease of DLA metallicity with increasing redshift 
(Prochaska et al. 2003; Rafelski et al. 2012) and a large in¬ 
trinsic dispersion of 0.5 dex out to z 5. Interestingly, 
this dispersion is not caused by observational uncertainties 
and does not appear to evolve with redshift. Neeleman et al. 
(2013) suggested that this effect could be partially explained 
by the existence of a fundamental plane in the redshift-mass- 
metallicity space, so that the dispersion in the metallicity- 
redshift relation is caused by the scatter in the masses of 
the dark matter (DM) halos hosting the DLAs (see also 
Ledoux et al. 2006; Mpller et al. 2013). There could be, how- 
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ever, additional sources of intrinsic scatter related to the en¬ 
vironment of the galaxy or the distribution of neutral clouds 
within the halo. 

DLAs have been extensively studied with hydrody¬ 
namic simulations (e.g. Nagamine, Springel & Hernquist 
2004; Pontzen et al. 2008; Fumagalli et al. 2011; Bird et al. 
2014), which have been able to reproduce the observed DLA 
abundance and metallicity range. While hydrodynamic sim¬ 
ulations are becoming extremely accurate and are now able 
to explain many of the properties of high-redshift galax¬ 
ies and their surroundings, their high computational cost 
prohibits studies of the large volumes of the Universe nec¬ 
essary for obtaining good statistics. Therefore the study 
of DLA metallicity dispersion with numerical simulations 
is extremely challenging. Notably, Cen (2012) simulated 
one overdense and one underdense region, whose properties 
bracketed many of the observed properties of DLAs. 

On the other hand, the semi-analytic approach, success¬ 
fully employed in the study of galaxy evolution, enables one 
to reproduce large galaxy populations and perform a com¬ 
prehensive parameter study at the cost of making several 
simplifying assumptions. Berry et al. (2014) studied DLA 
properties in the context of a detailed semi-analytic model 
of galaxy formation and were able to reproduce several key 
properties of DLAs, including the column density distribu¬ 
tion and the evolution of metallicity with redshift. However, 
these earlier models have so far just provided mean chemical 
evolution trends with z. 

In this work, we build a computationally efficient semi- 
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analytic merger tree model specifically designed to study 
the dispersion in the metallicity-redshift relation and the 
history of chemical enrichment of the Universe. We use the 
DM merger tree building algorithm from the GALFORM 
code (Cole et al. 2000) to construct regions with different 
histories of structure formation. We then apply the chemi¬ 
cal evolution model of Daigne et al. (2004, 2006) to explore 
the evolution of ISM metallicities in these different environ¬ 
ments. 

This paper is organized as follows. In Section 2 we de¬ 
scribe our model, which uses analytically computed DM halo 
merger trees and an accurate model of chemical evolution. 
In Section 3 we show that this model produces a dispersion 
in the observed DLA metallicities. We discuss our results in 
Section 4. 


2 MODEL 

Our goal is to calculate the dispersion in the mean mass 
fraction of collapsed structures, escape velocity and cos¬ 
mic star formation rate in different parts of the Universe. 
As a first step, we construct merger trees of DM halos us¬ 
ing the GALFORM algorithm (Cole et al. 2000). This al¬ 
gorithm relies on the extended Press-Schechter theory, in 
particular the conditional mass function (MF) of DM halos, 
and provides a list of mergers, and the redshifts at which 
they occurred, that led to the formation of a given halo. 
The modihcation by Parkinson, Cole & Helly (2008) ren¬ 
ders the algorithm consistent with the Sheth-Tormen MF 
(Sheth & Tormen 1999). Due to the shape of the matter 
power spectrum in cold DM models, the probability of a 
given halo with mass Mo to merge with another halo with 
mass M diverges as M —>■ 0, therefore we employ a cut¬ 
off mass Mres, which sets the mass resolution of our cal¬ 
culation. We use Mres = MminpO where M-min is the 
smallest halo mass that is able to form stars. In what fol¬ 
lows, we explore the range Mmin = 10^ — 10®Mq//i, and 
discuss this choice below. More details on the merger-tree 
building algorithm can be found in Cole et al. (2000) and 
Parkinson, Cole & Helly (2008). 

We take Vtot = 10® (Mpc/h)® as the total comoving 
volume of our calculation, having checked that for this vol¬ 
ume, the halo MF produced with the merger-tree algorithm 
converges to the Sheth-Tormen MF at all redshifts, and pop¬ 
ulate it with DM halos as follows. First we divide the total 
mass range M £ (M^es — 10 ^®Mq) into logarithmically equal 
bins of size A InM = 0.01. We then calculate the mean num¬ 
ber of halos in each bin at 2 = 0 using the Sheth-Tormen 
MF and draw the actual number N{M) from a Poisson dis¬ 
tribution with this mean. Finally, we sample N{M) halos 
with logarithmically equally distributed masses in the cor¬ 
responding bin. 

We build a merger tree for each halo at 2 = 0 and 
follow its evolution backwards in time up to Zf = 15, sav¬ 
ing the output in 50 redshift bins equally spaced in red- 
shift. In this manner, we obtain the distribution of DM ha¬ 
los in the comoving volume Vtot as a function of redshift. 
We further divide our sample into N — 1000 smaller regions 
AUi = Vtot/N = 10® (Mpc/h)® by assigning each merger 
tree to a random region i. This is equivalent to smooth¬ 
ing the density held on the scale of ro ~ 10 Mpc/h and 


assuming there is no correlation between neighbouring re¬ 
gions. This value was chosen because on scales larger than 
ro, galaxies are only weakly clustered (e.g. Norberg et al. 
2001; Zehavi et al. 2005), and we expect no signihcant dis¬ 
persion in the properties of different regions in the universe. 
On the other hand, for scales much smaller than ro there 
could be how of matter between different regions so that they 
cannot be treated as closed boxes. In other words, here we 
explore the effect of dispersion on the scale of massive clus¬ 
ters and deep voids; further contribution is expected from 
galactic scales, as we briehy discuss below. For comparison, 
in the simulations of Cen (2012) the overdense and under- 
dense regions had a size of roughly 20 and 30 Mpc/h on the 
side, respectively (see also Rollinde et al. 2013). 

We then calculate the mean mass fraction in collapsed 
structures fcoii, the mean escape velocity and the mean star 
formation rate (SFR) in each region AVt as explained be¬ 
low. While the mean MF and fcoii in our total volume Vtot 
coincide with the expectations from the Sheth-Thormen MF 
at all redshifts, their values differ in each AU; region. Those 
regions that host a group or a cluster at low redshift have a 
higher concentration of structures already present at higher 
redshift, whereas present day voids are relatively empty at 
early times. This small-scale inhomogeneity creates the ob¬ 
served dispersion in the different observables. 

The mean mass fraction in collapsed structures fcoii is 
calculated by summing the masses of all the halos above 
Mmin in the given region and dividing by AVi and the mean 
matter density pm'- 


fcoii,i 


AVi Pm 


( 1 ) 


The dispersion in fcoii among the different regions grows 
with redshift for all values of Mmin, reaching 0.35 dex for 
Mmin = IO^Mq/H at 2 = 10 and 0.25 dex at 2 = 5. Fur¬ 
thermore, the dispersion grows with Mmin, since if only large 
structures are allowed to form, some regions remain practi¬ 
cally devoid of structures. In the case of Mmin = 10^ M^/h 
and Mmin = 10^Mq//i the dispersion is about 0.25 dex and 
0.18 dex at 2 = 10, respectively. 

The mean square of the escape velocity Vcsc{z) within 
each region is calculated as the mass-weighted average of 
2GM/R for each halo mass M, where R is the virial ra¬ 
dius enclosing a volume with mean density 200 times the 
critical density of the Universe. Contrary to the case of 
the collapse fraction, the dispersion in Vcsc is slightly larger 
for smaller values of Mmin, being on average 0.2 dex for 
Mmin = Vp Mq/H. The reason for this is that models with 
low Mmin exhibit a larger range of masses in each region, 
so that while the total mass in two representative regions 
might be the same, their mean escape velocities could differ 
substantially. In the case of high Mmin, however, the mass 
spectrum of different regions is much more uniform and the 
main difference is in the total mass of collapsed structures, 
reflected in fcoii- Also contrary to the case of fcoii the disper¬ 
sion in Vcac grows with decreasing redshift. The mean escape 
velocity is used in our chemical evolution code to calculate 
the outflow rate as explained below. 

Our next step is to calculate the star formation rate 
(SFR) in each halo and the mean cosmic SFR in each region 
AVi- We use the results of Behroozi, Wechsler & Conroy 
(2013) (hereafter B13) to calculate the SFR as a function of 
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halo mass M and redshift z in the range 0 < 2 < 8 and 9 < 
logjQ(M/M 0 ) < 15. We extrapolate these results to lower 
masses and higher redshifts with SFR oc 
where j3 = 2.1 was obtained from a fit to the results of 
B13 in the range 5 < 2 < 8 and 9 < logj^g(M/M 0 ) < 11 
and a = —0.1 was obtained from the observed cosmic SFR 
at 2 ~ 10 (Bouwens et al. 2014; Oesch et al. 2014). We note 
that this extrapolation is purely phenomenological and IMF- 
independent. We calculate the SFR for each halo and then 
compute the cosmic SFR density ipi{t) as a function of red- 
shift in each region AVi: 


= AT/. - 


( 2 ) 


We note that by using the fit from B13 (which refers to the 
universal SFR density) in all of our calculations, in particu¬ 
lar the small-scale regions, we implicitly assume that the star 
formation efficiency depends only on the DM halo mass and 
redshift. This approximation is justified by the fact that we 
compute mean values over sufficiently large regions, which 
we treat as closed boxes. The resulting dispersion in the 
SFR between different regions is roughly 0.4 dex, dropping 
significantly below 2 ~ 2.5. 

Having obtained the rate of growth of structure and the 
SFR density in each region we proceed to the calculation of 
the chemical enrichment of the interstellar matter (ISM). We 
use the chemical evolution model developed in Daigne et al. 
(2004, 2006) and Rollinde et al. (2009) which follows the 
exchange of mass between the gas within and outside of 
collapsed structures, the SFR at each redshift and the rate 
of metal production in stars. 

The initial gas content of galaxies is taken to be equal 
to the cosmic mean fbaryon = The model then fol¬ 

lows two gas reservoirs, the intergalactic matter (IGM) and 
the ISM. Matter flows from IGM to ISM as the galaxies 
form, where baryons are assumed to follow DM, so that 
the mean baryon accretion rate in each region is given by 
ab{t) = ^Ibdfcoii/dt and fcoii is calculated as in eq. (1). Once 
inside galaxies, baryons form stars according to the rate ip{t) 
calculated above (see eq. (2)). We assume a Salpeter initial 
mass function (IMF) 4>(m) for m™/ < m < rUsup with 
rriinf = O.IM 0 and rrisup = IOOM 0 . 

Baryons can flow from structures back to the IGM due 
to galactic winds or feedback from SNe or AGN. In this work 
we consider outflows powered by stellar explosions with the 
rate given by: 


o(t) = f dm^{m)^{t-T{m))Ekiri{rn) (3) 

vLae J 

where r(m) is the lifetime of a star with mass m, Ekin 
is the kinetic energy released when this star dies, Vesc is 
the mean escape velocity from structures and e = 0.0006 
is the fraction of kinetic energy that powers the outflow. 
This value was chosen as in Daigne et al. (2006) through 
comparison of the model predictions for the IGM compo¬ 
nent and the observed oxygen and carbon abundances in the 
Lya forest (e.g. Simcoe, Sargent & Rauch 2004; Songaila 
2005; Aguirre et al. 2008). We verified that small deviations 
from this fiducial value do not have a significant effect on 
the resulting ISM metallicity. The contribution of this term 
strongly depends on Mmin, since small halos have smaller 
escape velocities and as a result are more easily disrupted. 


0.5 
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Figure 1. Log of the metallicity abundance relative to the so¬ 
lar value for 100 regions, each with a volume of AVi = 10® 
Mpc®/h® (thin grey lines) and their mean (thick blue line) for 
Mmin = 10 ^Mq/H. The dispersion among the realizations and 
their mean is represented by the blue points with error bars. The 
mean dispersion in the redshift range 2 = 2 — 5 is 0.25 dex. Red 
dashed lines represent the estimated upper and lower limits when 
the mass-metallicity relation within each region is also considered 
(see text). Black points represent data from Rafelski et al. (2012). 


Consequently, models with lower Mmin will be able to retain 
fewer baryons inside the galaxies and will have lower ISM 
metallicities, as shown below. 

To sum up, the evolution of the baryonic mass in the 
IGM and ISM is given by dMioM/dt = —ab{t) + o{t) and 
dMisM/dt = [—'ip{t) + e(f)] -I- [ab{t) — o{t)] where e{t) is the 
rate at which stellar mass is returned to the ISM by mass 
loss or stellar deaths. These equations are solved for each 
region separately, so that no matter flow is allowed between 
regions. This approximation is justified by our choice of the 
region size. 

The chemical evolution of the IGM and ISM is cal¬ 
culated as described in Daigne et al. (2004). In particular, 
we do not employ the instantaneous recycling approxima¬ 
tion but calculate the rate at which gas is returned to the 
ISM including the effect of stellar lifetimes and computing 
stellar yields for each element and for different stellar mass 
ranges. Further details on the chemical evolution model can 
be found in Daigne et al. (2004, 2006) and Vangioni et al. 
(2015). 

In the next Section we show results for the metallicity of 
the ISM relative to the solar value: [M/H] = logjQ(M/H’) — 
logjQ(M ///)0 and compare it to DLA data. 


3 RESULTS 

We can now test our model against observational data, using 
DLAs as probes of the ISM at high redshift. For our fiducial 
model we take Mmin = llfMQ/h. We used 100 regions, 
each with AVi = 10® (Mpc/h)®, where both the mean and 
the variance (calculated in redshift bins of A 2 = 0.5) over 
the whole ensemble converge. 

In Figure 1, we show the evolution of the metal abun¬ 
dance relative to the solar value for our fiducial model for 
100 regions (grey lines) and compare them to observations of 
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DLAs from Rafelski et al. (2012). The thick blue line shows 
the mean of all the regions and the error bars represent their 
variance. It can be seen that the dispersion predicted by the 
model (0.25 dex in the redshift range of the observations) is 
significant, but smaller than the observed dispersion in the 
metallicity-redshift relation (0.5 dex). We can thus conclude 
that the dispersion in the fraction of collapsed structures, 
escape velocity and SFR in each region AVi contributes 
to the dispersion in the metallicity-redshift relation of the 
DLAs, but is not its only source. We note that the disper¬ 
sion predicted by our model is a lower limit on the actual 
value, since the observed metallicity of a given DLA system 
might also depend on the mass of the host DM halo. Indeed, 
Neeleman et al. (2013) showed that by applying a correction 
following from the existence of a global mass-metallicity- 
redshift plane, the scatter in the metallicity-redshift relation 
reduces from 0.5 dex to 0.38 dex. In addition, the potential 
metallicity gradients of DLAs can also contribute to the dis¬ 
persion (Chen et al. 2005; Christensen et al. 2014). Keeping 
in mind that the observational uncertainty for each mea¬ 
surement is about 0.12 dex it is tempting to suggest that 
the combination of a mass-metallicity-redshift plane and the 
effect explored here, which stems mainly from the different 
formation epochs of structures in over- and underdense re¬ 
gions, might explain the observed dispersion. We stress that 
although we expect a correlation between these two effects, 
they are somewhat distinct. Indeed, the dispersion in the 
mean log Mhaio between regions is below 0.01 dex for z < 5. 

While the modeling of the full effect is beyond the scope 
of the present paper, we tried to estimate it using the ob¬ 
served mass-metallicity relation from Maiolino et al. (2008) 
and the dispersion in stellar masses within each region using 
our merger tree, assuming for simplicity a constant stellar 
mass to halo mass ratio. Since the mass-metallicity relation 
is redshift related, and the parameters for high redshifts are 
somewhat uncertain (see also Zahid et al. 2013), we used the 
relation for 2 = 3.5 for the whole redshift range shown on 
Figure 1. Adding the resulting dispersion to the metallicities 
calculated with our model (grey lines) produces the upper 
and lower limits shown as the red dashed lines on Figure 
1. It can be seen that this very rough estimate can indeed 
explain the whole obserbed range of DLA metallicities, al¬ 
though we caution that a more self-consistent treatment is 
needed, which we leave to future work. 

The choice of the minimal halo mass that is able to 
form stars affects both the normalization and the disper¬ 
sion of the predicted metallicity, as can be seen in Fig¬ 
ure 2 where we show model results for Mmin = 10^, 10® 
and IO^Mq/H (from bottom to top). In the case of high 
Mmin galaxies start forming later, but they are on aver¬ 
age more massive and have higher escape velocities. As 
a consequence, the outflows are weaker and a larger frac¬ 
tion of metals produced in stars stay in the ISM as com¬ 
pared to the case of low Mmin- It is interesting to note 
that the explored range Mmin = 10^ — ItEM^/h matches 
the mean of the observations. We note that the numeri¬ 
cal simulations of Nagamine, Springel & Hernquist (2004) 
show that the halo mass scale below which no DLAs 
exist is slightly above 10®M© at 2 = 3 — 4, whereas 
Pontzen et al. (2008) find that the occurrence rate of 
DLAs drops sharply below M 10^ Mq (although see 
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Figure 2. Log of the metallicity abundance relative to the solar 
value, averaged over 100 regions. Error bars show the dispersion 
for each model with Mmin = 10^, 10® and ItP Mq/H (from bot¬ 
tom to top). Black points are measurements from Rafelski et al. 
(2012). The inset shows the dispersion (in dex) for Mmin = 
10^, 10® and IO^Mq/Zi (from top to bottom). 

Webster, Bland-Hawthorn & Sutherland 2015). Our results 
are also consistent with these limits. 

The inset in Figure 2 shows the dispersion in metallic¬ 
ity among 100 regions as a function of redshift. It can be 
seen that the dispersion is larger for smaller Mmin- In fact, 
there are two competing effects: in models with larger Mmin 
the dispersion in fcoii is larger, but the dispersion in is 
smaller. While fcoii determines the overall number of galax¬ 
ies in each region, influences the outflow, and is partic¬ 
ularly important in low-mass halos which could thereby lose 
a substantial fraction of their newly acquired metals. 


4 DISCUSSION 

In this paper, we have used the chemical evolution model of 
Daigne et al. (2004, 2006), supplemented by a merger-tree 
description of structure formation and an observation-based 
estimate of the SFR, to study the metal enrichment process 
of the ISM. Our models include the full range of galaxy 
formation physics required to study the chemical evolution 
of DLAs, including gas accretion, star formation rates, star 
formation histories, gas outflows and mergers. We reached 
two important conclusions: 

(i) The dispersion in structure formation rates in dif¬ 
ferent parts of the Universe and the accompanying variance 
in the stellar content of protogalaxies produces a dispersion 
in the metallicity-redshift relation of DLAs of roughly 0.25 
dex. Furthermore, this dispersion depends on the assumed 
minimal mass of halos that are able to form stars Mmin 
through the dispersion in the structure formation histories: 
models with smaller Mmin produce larger dispersion. 

(ii) The range Mmin = 10^ — lOM^/h provides a good 
description of the data with little variance between these 
models. This mass range corresponds to halo masses with 
virial temperatures above T > 10"^ for which neutral gas 
cooling is inefficient. 

The dispersion resulting from our model is smaller than 
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the observed value of ~ 0.5 dex, but is clearly an impor¬ 
tant ingredient in the overall effect. Another potentially im¬ 
portant source of intrinsic scatter stems from the difference 
in the masses of DM halos hosting the DLAs, for exam¬ 
ple through the global mass-metallicity-redshift relation dis¬ 
covered by Neeleman et al. (2013). Specifically, even within 
each over- or underdense region, there is a certain spectrum 
of halo masses with different metal enrichment histories. Ac¬ 
cording to our very rough esimate, the combination of these 
two effects can explain the observed dispersion, but a more 
rigorous treatment is needed, which takes into account their 
possible correlations. Another important effect is related to 
the structure of the galactic halos giving rise to the DLAs, 
which may also contribute to the observed metallicity dis¬ 
persion. 

Although the conclusions outlined above are relatively 
robust, our model involves several approximations which can 
be readily improved. One such assumption is the SFR as 
a function of halo mass and redshift, in particular our ex¬ 
trapolation of the results of B13 to high redshifts and low 
masses, down to M = IO^Mq. It would be beneficial to 
explore more accurate and physically motivated extrapola¬ 
tion schemes, for example using full semi-analytic models of 
galaxy evolution. 

Galactic outflows obviously play an important role in 
galaxy evolution, and in particular in the process of metal 
enrichment of the ISM. We have partly addressed this is¬ 
sue by considering different values of Mmin which resulted 
in different mean escape velocities and outflow rates. This 
effect is crucial in setting the correct metallicity normal¬ 
ization, and we have showed that models with lower Mmin 
have lower mean escape velocities and retain less metals in 
the ISM. Clearly, a more comprehensive study of the effects 
of different feedback mechanisms on the normalization and 
dispersion in DLA metallicities is warranted, and we leave 
this to future work. 

Finally, we plan to extend the model presented here to 
a full description of chemical enrichment on a halo-to-halo 
basis. This will allow us to address the issue of a global 
mass-metallicity-redshift relation of DLAs, as well as the 
effect of the impact parameter of quasar sightlines, thereby 
accounting for the different sources of metallicity dispersion. 
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